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Abstract 



We model the evolution of the concentration field of macromole- 
cules in a symmetric field-flow fractionation (FFF) channel by a one- 
dimensional advection-diffusion equation. The coefficients are pre- 
cisely determined from the fluid dynamics. This model gives quanti- 
tative predictions of the time of elution of the molecules and the width 
in time of the concentration pulse. The model is rigorously supported 
by centre manifold theory. Errors of the derived model are quantified 
for improved predictions if necessary. The advection-diffusion equa- 
tion is used to find that the optimal condition in a symmetric FFF 
for the separation of two species of molecules with similar diffusivities 
involves a high rate of cross-flow. 
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1 Introduction 

Consider the transport of some contaminant molecules in the fluid flow of 
a symmetric field-flow fractionation (FFF) channel as analysed by Giddings 
and others |5], 13, e.g.] and sketched in Figure^. The two horizontal parallel 



plates above and below the channel are not permeable to the contaminant 
molecules but allow for the cross-flow of fluid. This cross-flow distributes 
the contaminant preferentially to the lower side of the channel as shown in 
Figure [| It is this cross-flow and asymmetric distribution of contaminant 
concentration c(x, y, t) that creates a differential advection of different molec- 
ular species and renders the problem interesting. 

Using techniques based upon centre manifold theory fp70| , from the con- 



tinuum equations (Section we deduce that a model for the contaminant 
distribution in the channel is the advection-diffusion equation 

— = -U— + D— (1) 

dt dx dx 2 ' 

where t denotes time, x measures distance downstream along the channel, and 
C(x,t) = c(x,0,t) is the concentration of the contaminant measured along 
the lower plate (the so called accumulating wall). We derive expressions for 
the effective advection velocity U as it predominantly determines the time of 
efflux of the contaminant out across the end of the channel, and the effective 
diffusivity D as it determines how wide the contaminant spreads by the time 
it reaches the end of the channel: in a useful parameter regime (Section |) 

^ 6uk p 72u 2 k 3 
v b ' v^b 2 ' 

where k is the molecular diffusivity, u is the mean along-channel velocity, b is 
the channel height, and vq is the cross-flow velocity. The term -D§^ models 
the so called "zone broadening effects" discussed by Litzen and others J7L [T3 . 
We also quantify the two sources of errors in the model by 

• estimating the time it takes for initial transients to die out and the 
model to become valid (Section Rf); 
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Figure 1: Side view of symmetric field-flow fractionation (FFF) channel. 



• determining higher-order corrections to the advection-diffusion model 
(Section §. 

This model and its errors may be rigorously justified as discussed in other 
applications of centre manifold theory to shear dispersion by Mercer, Roberts 
and Watt B g g 0, . 




Field-flow fractionation channels are used to separate species of contam- 
inant molecules with different diffusivities. In Section [5| we use model ([!]) to 
identify that FFF separates molecular species most efficiently for relatively 
high cross-flow: up to 

(3) 

Consequently, in describing the governing equations in Section || we introduce 
a non-dimensionalisation appropriate for such high cross-flow rates. 

Further research in field-flow fractionation could model the dynamics of 
contaminant molecules in tubular channels ||14|| , trapezoidal channels f7], or 
in asymmetric FFF channels 0, as well as the dynamics of non- neutrally 
buoyant particles |L3 . 



2 Governing equations for symmetric FFF 

Consider a symmetric FFF channel as discussed by Giddings and others 
113], T7] and depicted schematically in Figure |T|. The dynamics takes place 
between two flat plates located at y = and y = b. The fluid flow between 
the plates is driven predominantly by a pressure gradient p x parallel to the 
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Table 1: Typical set of physical parameters for FFF and the consequent 
parameters (in the second part) appearing in the analysis. The data is for 
the Cow Pea Mosaic Virus g p464] in the FFF channel of |TJ] . 



Parameter 


Value 


Channel width 


b 


0.05 cm 


Kinematic viscosity 


V 


0.01cm 2 /s 


Mean longitudinal velocity 


u 


0.1 cm/s 


Cross-flow velocity 




5xl0~ 4 cm/s 


Molecular diffusivity 


K 


2 x 10" 7 cm 2 /s 


Boundary layer (BL) thickness 


V 


4 x 10~ 4 cm 


Cross-BL diffusion time 


T 


0.8 s 


Longitudinal velocity in the BL 


U 


5 x 10~ 3 cm/s 


Downstream advection distance 




4 x 10~ 3 cm 


Prandtl number 


a 


5 x 10 4 


Cross-channel Peclet number 


V 


125 


Downstream Peclet number 


u 


2.5 x 10 4 


Velocity ratio 




0.1 



plates. Being that of a Newtonian fluid with kinematic viscosity v and density 
p, the velocity field is essentially that of parabolic Poiseuille flow except that 
there is a cross- flow, of velocity — Vq, from the upper plate to the lower (if v 
is positive) . The plates are permeable to the fluid in order for this cross-flow 
to occur; but they are impermeable to the contaminant molecules. Within 
the fluid the contaminant, of concentration c(x,y,t), is advected by the flow 
and diffuses with coefficient k. In this section we non-dimensionalise the 
governing differential equations, and also deduce the advecting fluid velocity 
field and confirm that it is nearly parabolic. 

For order of magnitude estimates of quantities we use the geometry of 
Wahlund Sz Giddings [ 14| : the channel width is b pa 0.05 cm; the density of 
the fluid, water, is p = 1 gm/cm 3 ; and the kinematic viscosity v pa 0.01 cm 2 /s. 
The fluid moves so that on average it takes about 5-15 minutes to traverse 
about 50 cm so a typical fluid velocity is u ~ 0.1 cm/s and the driving pressure 
gradient must be roughly p x pa 5gm/cm 2 /s 2 . The cross-flow is driven at 
rates vq of order 5 x 10 _4 cm/s. When the contaminant molecules are the 
Cow Pea Mosaic Virus ||, p464], this configuration gives parameters as listed 
in Table |I|. We base our analysis on this set being typical of parameters of 
interest. 

The equations governing the fluid motion are the Navier-Stokes and con- 
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tinuity equations 

-r^ + q- Vq= — Vp + z/V 2 q, (4) 
ot p 

V ■ q = (5) 

for the incompressible velocity field q = ui + v j and for the pressure p. The 
contaminant evolves according to the advection-diffusion equation 

dc 

— + q-Vc = KV 2 c (6) 
ot 

for the concentration field c. Herein we assume the molecules are neutrally 



buoyant, though sedimentation effects JT3|] could be included in further work 
by modifying this equation. Note that although we are concerned with 
the dynamics of the concentration field c, we only seek the steady and x- 
independent fluid flow governed by the Navier-Stokes and continuity equa- 
tions. The boundary conditions on the plates are those of no longitudinal 
flow, 

u = , v = —vq , on y = and y = b, (7) 
and no flux of the contaminant through the plates, 

dc 

VqC + k— — = , on y = and y = b . (8) 
dy 

The above equations fully specify the dynamics of the fluid and the contam- 
inant molecules in the channel. 

The non-dimensionalisation we adopt is chosen to reflect the fact that for 
the regime of most effective separation of species (see Section ||) the contam- 
inant is concentrated near the lower plate due to the cross-flow. Introduce 
the following non-dimensional variables denoted by stars: 

x =-, y = - , t = -, u = — , v = — , P = — , (9) 

where rj = k/vq is the characteristic thickness of the distribution of contam- 
inant in a boundary layer near the lower plate, r = tj/vq = if j k = k/vq is 
the cross-boundary layer advection (t]/v ) or equivalently the cross-boundary 
layer diffusion {rf / ' k) time, 

1 dp bin Qu 

2 ox pv V 

is the characteristic downstream velocity in the boundary layer, u = — ^ 
is the mean speed of the Poiseuille flow in absence of the cross-flow, £ = u r 
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is the downstream advection distance for the material in the boundary layer 
in a cross-boundary layer diffusion time, and where 



— and V = 

K K 



(n; 



are Prandtl and cross-channel Peclet numbers, respectively. Typical values 
of all these quantities are recorded in Table |l|. In essence this scaling is that 
of the distribution of contaminant molecules which typically are swept to be 
near the lower plate with the upper plate "far away" at y = V. Substitute 
these scalings into the equations and omit the distinguishing stars hereafter. 

The steady fluid flow is straightforward to determine. The ^-momentum 
equation determines that v — — 1 everywhere. The x-momentum equation 
for the steady velocity field u(y) becomes 



V 

2 



1 du d 2 u 
o dy dy 2 



+ 1 = 



(12) 



with boundary conditions u(0) 
velocity component is 



u(V) = 0. The exact solution for this 



u(y) 



V- 



-y/v 



-V/<T 



y 



' y 2 \^y V y+y 2 \+ n v 2 

y ~ V ) + o U"2 + 3V + °b 



(13) 
(14) 



Observe, as used in earlier analyses [jM, e.g.], the downstream advection is 
nearly parabolic; because the Prandtl number a is so large the correction of 
O (V/cr) is usually negligible. 

The dynamics of the contaminant remains nontrivial. Under our nondi- 
mensionalisation the advection-diffusion equation becomes 



dc dc dc 
dt dx dy 



d 2 c ^ 2 d 2 c 



dy 



dx 2 



(15) 



where 



K 



and U 



ub 



(16) 



uq 6 u QU k 
are respectively the velocity ratio and the downstream Peclet number based 
on the mean longitudinal speed, see Table ffl. The non-dimensional boundary 
conditions for the contaminant are 

dc 



c + 



dy 



at y = and y = V . 



(17) 
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We analyse the dynamics described by this non-dimensional equation in 
this paper. The main non-dimensional parameter V, appearing as the non- 
dimensional width of the channel, is typically large, O(10 2 ), as we expect 
cross-flow advection to keep the contaminant close to the bottom plate. 



3 The dynamics approach a centre manifold 

We justify the basis of model ([!]) using centre manifold theory JTJ as adapted 



10, U to the long thin geometry of the FFF channel. Under the action of 
the cross-flow balanced by diffusion the contaminant distribution across the 
channel relaxes quickly to an exponential distribution, c = Cexp(— y). The 
shear velocity, different at different y, will smear this contaminant cloud out 
along the channel, while cross-flow and diffusion continue to act to push the 
cross-channel distribution towards the exponential distribution. The net ef- 
fect is that the cloud has a concentration that is slowly varying along the 
channel and is approximately exponential across it. Thus, after the quick 
decay of cross-stream transients, we justify the relatively slow long-term evo- 
lution of a contaminant cloud for which x derivatives of C, d n C/dx n , are 
small. 

An initial "linear" picture of the dynamics is established by assuming 
that there are no downstream variations. When downstream gradients are 
ignored, the relaxation across the channel of the contaminant obeys the dy- 
namics 

dc dc d 2 c dc . _ /10 . 

^7~^ _ = ^> s - t - c+^- = 0ony = 0and?/ = V. (18) 
at ay ay 1 ay 

The neutral solution already mentioned is the exponential Cq = Cexp(-y). 
The other solutions, all decaying, are 



a 



■ ■\~-\ l t) ■ (19) 



for n = 1,2,..., where C n are constant coefficients determined by the initial 
condition such that 



c(y,0) = ^c n (y,0) (20) 



n=0 



and the decay rate is 



1 n 2 ix 2 . . 

n = ~4 y2~ • (21) 

The slowest rate of decay to the centre manifold will be due to the 
n = 1 mode, although in many cases the second term in (|21|) is negligi- 
bly small as V is of order 10 2 . As an example, assume an initially uniform 
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distribution of contaminant across the channel, then from ( p0|) expect that 
C\ oc exp(V/2). The relaxation process is then dominated by the exponen- 
tial decay of exp(V/2 + Ait) which effectively leads to a relaxation time of 
roughly t re i = — V/(2Ai). In dimensional form this cross-channel relaxation 
time 

t rel « 2Vr = 200 s , 

which agrees with the experimental observations of several minutes for macro- 
molecules given, for example, in |b|, Eqn (30)]. Thus expect the decay to a 
low-dimensional centre manifold to occur on this time scale. 

The presence of downstream x-variations perturbs the contaminant pulse 
and results in its non-trivial long-time evolution. Centre manifold theory pro- 
vides a powerful rationale for modelling such evolution where the long-term 
behaviour is separated from rapidly decaying transients. This was recognised 
by Coullet & Spiegel |§] and Carr & Muncaster [§, U; see the draft review 
by Roberts JTT] for an extensive discussion. The application of the theory to 
dispersion in channels and pipes has been developed by Roberts, Mercer and 



Watt pi RL P, 16, 15|. Using the same techniques here, we seek a solution 



to the governing equations in the form 

3C 

c = h(C, y) such that — = g(C) . (22) 

Here the function h, C exp(-y) to leading approximation, describes the de- 
tails of the contaminant field throughout space and time in terms of the 
concentration C of contaminant at the lower plate. Such a solution forms 
a model of the dynamics for two reasons. First, the low-dimensional set of 
states described by h(C) are exponentially attractive because of the action of 
cross-stream advection and diffusion as seen above. Secondly, the associated 
function g models the effective advection and diffusion of the contaminant in 
the horizontal by describing the evolution of C. 

We find approximations to these functions by assuming that the con- 
centration field is slowly varying in the horizontal, that is, d/dx is a small 
operator. Rigorously, one would expand in the downstream wavenumber as 



introduced by Roberts fllOl . Formally we express h and g in the following 
asymptotic series 

9~Y.9n-r— and h ~ £ h n (y) — - , (23) 

n=l UJj n=0 UX 

where for example h = exp(-y) is the leading order approximation to the 
contaminant field, —g\ = U is the effective advection velocity, and g<i = D is 
an effective horizontal diffusion coefficient. The advection-diffusion model ([!]) 
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is obtained from just the first two terms in the expansion for g. In dispersion 
problems, the asymptotic series in (|23|) typically converge in a sense discussed 
by Mercer, Roberts and Watt §, g [LjJ. 

To find the asymptotic expansions (|23D we implement an iterative al- 
gorithm (see [|l^]) in computer algebra (see Appendix |A|). The results are 
assured to be accurate by the approximation theorem of centre manifolds. As- 
sume that some approximate solution of the contaminant advect ion-diffusion 
equation ( |T5| ) with boundary conditions ( |T7| ) is found in the centre mani- 
fold form (|22|) ; for example, the iteration is initiated with the approximation 
c = C exp(-y) and g = 0. We wish to refine such an approximation by 
finding a correction h! to the shape of the centre manifold and a correction 
g' to the evolution thereon. As established by Roberts |T2J the corrections 
are found by solving 

d 2 ti dti 

w + t = R+3 ' eM - vh (24) 

where R is the residual of ([H|), with boundary conditions 
dh! 

ti + — = at y = and y = V, and U = at « = 0. (25) 
dy 

This last boundary condition reflects that we seek a solution parameterized 
by the concentration at the lower plate: C(x,t) = c\ y=0 . The correction to 
the evolution g' is chosen to satisfy the solvability condition 

[ V R + g'ex V (-y)dy = (26) 
Jo 

in order to satisfy boundary conditions (EH). Then the differential equa- 



tion (|24 ) is solved to find h! . The iterations continue until the desired terms 
are found in the asymptotic approximation to the centre manifold (p3|). Com- 
puter algebra, such as the program listed in Appendix |A|, easily performs all 
the algebraic details. 



4 The detailed centre manifold model 

Since all the algebraic machinations are handled by the computer algebra 
of Appendix 0, here we just record and discuss the results. General results 
simplify considerably in the typical case of large V when the contaminant is 
held near the lower plate. Then higher order corrections are readily found. 
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From the computer algebra results, the concentration field (the centre 
manifold) is to low order 



Ce- y 
_ dC_ 
dx 

1 dC e~ y y 2 
~ 12o~~dx~ V 



i-I 

a 



2me 



-v 



y 



e- y 2me- v (y + 1)- — + + ^ 



V 3V 



{V-yf + O 



' d 2 C 
dx 2 ' a 2 



(27) 



where the evolution of the contaminant concentration along the bottom plate 
is described to leading order by 



d_C_ _ dC_ 

dt dx 



i-I 

a 



2\ V 

l - 2m+ v)-lTo 



o(— — 

\ dx 2 ' a 2 



(2* 



(29) 



where 

m= ( 1 - 6 ) ~{ l/V + 1/2 asV-0 ' 

The order of error notation O (a, /3) is used to denote errors O (a) + O {(5). 

Since the typical cross-channel Peclet number V is of order 10 2 we take 
m = 1 in presenting further detailed results (for completeness we present 
results for weak and moderate cross- flows in Appendix [B|). The dominant 
error in this approximation is O ( e_V ) an d so expect it to be acceptable for 
V greater than about 6. Then (1271) simplifies to 



Ce^ + ^y 2 
ox 



+ 



' d 2 C 
dx 2 ' a 2 



i-I 

a 

-VI 



3 + y_ 
3V 



12V a 



(30) 



This shows the predominantly exponential distribution of the contaminant 
as advection towards the lower plate by the cross-flow is counter balanced 
by diffusion. The exponential distribution is modified by the interaction of 
the shear flow and the along-channel spatial gradients of the contaminant 
as given by the second term in (^0J) and shown in Figure 0. The above 
expressions give the details of the concentration field parameterized by its 
value C(x,t) = c\ y= o at the lower plate. 

The associated advection-diffusion equation is (|l|) with coefficients 
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Figure 2: Velocity field and an instantaneous concentration field near the 
accumulating lower plate (a) when the concentration along the wall is given 
by the Gaussian (b). Fields correspond to the parameters given in Table |I|. 
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2 .-, 20 , 56 \ 2(V - 8) ( 12 , 42^ 




(32) 

giving the effective advection speed and dispersion coefficient. The crudest 
approximation, but useful over a reasonable parameter regime, is that U ~ 1 
and D « 2 leading to the dimensional expressions given in the Introduction. 

Running the computer algebra program to higher order in spatial gradi- 
ents we find that the dynamics of the dispersion is governed by the extended 
evolution equation 

D ^£_ E & c F dA ° o(—\ (33) 

dt dx dx 2 dx 3 dx 4 \ dx 5 I ' 

where the coefficients of the third and fourth order derivatives are: 

/ 102 744 1936 \ , , 

E = - 4 ( 5 " — + ) < 34 > 

2 / „ 2220 13408 31856\ m 

__( 5 V-170 + — -— + W ) + C 

725 9480 58292 142168 

16 / , „ 44145 464560 2549556 5856960\ 
+ _(44V-2175 + ^ + —) 

(35) 

The d 3 term in (^) with coefficient governs the skewness of the pre- 
dictions of the model by modifying the effective advection speed of various 
spatial modes. The d% term with coefficient F affects the decay of the spatial 
modes. Note that F is positive for at least large enough V and a — a fourth 
order model may thus be unstable for short enough spatial modes: approxi- 
mately the fourth order model is unstable for along channel non-dimensional 
wavenumbers \k\ > 1/(4 \/IT). Thus although the third-order term may be 
used to improve predictions of the advection-diffusion model, the fourth-order 
term should be limited to helping estimate errors in the predictions. 





5 Species separate best at high cross-flow 

The aim of field-flow fractionation is to separate as far as possible two or 
more different species of contaminant molecules. Different contaminants are 
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characterised by different diffusivities, k 3 - say. A contaminant with lower 
diffusivity will be pushed closer to the lower plate by the cross-flow. Conse- 
quently, its effective advection speed along the channel will be lower. Thus 
one collects a contaminant with higher diffusivity at the exit before a contam- 
inant with lower diffusivity. Here we identify the operating regime when the 
separation is most effective between two species of nearly the same diffusivity. 

Consider the advection-diffusion predicted by model (|l|) with different 
species identified by the subscript j. In the non-dimensional analysis this 
leads to different characteristic scales: 



1l c 

v 



UojTj 



(36) 



Thus the advection-diffusion model ([J) for the jth species has dimensional 
coefficients 

s i 



Uj = u 0j U(Vj) , Dj = ^-D(V, 



(37) 



from the leading term in each of (Bjj)— (^3) upon neglecting terms of order 
Vj/cTj. In a channel of fixed length L the approximate times of efflux are 



U- 



L 



L Vi 



u 03 U{V 3 ) 6uU(Vj) 



(3* 



Then the time interval between the moments when the two contaminant 
pulses with diffusivities k± — k — Ak/2 and = K + Ak/2, injected simul- 
taneously at the beginning of the channel, exit the channel is 



AT 



dT 
Ok 



PL \U(V) — VU'(V) \ An 



6 u 



U 2 {V) 



K 



(39) 



The width of the contaminant pulse at the time of efflux is proportional to 
V DT, and hence the time taken for a contaminant pulse to pass the end is 
proportional to 5 where 



DT t 2 D(V) 1L b 2 D(V) 



L 



(40) 



U 2 ~ {, U 3 (V) 6ukVU 3 (V)' 

To maximise separation of two species with close values of diffusion coeffi- 
cients Kj we need to maximise the difference in the time of efflux relative to 
the width in time of the pulses at the efflux. Thus for a given small change 
in diffusivity, An > 0, we wish to maximise 



AT 


dTdV 


An 


V 


dT 


An 


5 




~T 


~ 5 


dV 


n 



L V 3 / 2 \U(V) - VU'(V)\ Ak 



QUA 



D{V)U{V) 



(41) 
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Expect the existence of an optimum cross-flow from the following physical 
arguments. Increasing the cross-flow significantly increases the difference 
in efflux times. On other hand, an extremely strong cross-flow would keep 
both contaminants close to the plate in a very slow flow for long enough so 
that longitudinal molecular diffusion becomes significant. Thus resolution 
will decrease for an excessively strong cross-flow. The optimal separation of 
species with given diffusivities in a channel of fixed geometry with a fixed 
fluid flux through it is accomplished when 

= v^\u(v)-vu'(v)\ 

'D(V)U(V) 

«^(V-4) (42) 
'(V - 2) (V 6 + 36W 2 (2V 2 - 20V + 56)) 

is maximised. From dR/dV = we obtain 

V 6 (V 2 - 12V + 16) = 72U 2 (3V 4 + 52V 3 - 336V 2 + 912V - 896) (43) 

with a solution for optimal V of 

^^+f-f^ +0 (I). 



The leading term of this optimum gives the optimum vq mentioned in the 
Introduction. As seen from Figure |3|, for the parameter values listed in 
Table |l]this optimum occurs at Vo ~ 613 which corresponds to the relatively 
high cross-flow velocity Vq ~ 2.5 x 10~ 3 cm/s. Then the optimal regime of 
two species separation gives 




AT 6 1 / 8 /3LA« 1/4 



24 384 




(45) 



For the geometry of the channel considered in [II] and parameters given in 
Table [l] the maximum resolution is thus 

AT A/t 

— « 430 , 46 

K 

where 8 ~ 1 min. For the regime considered the time necessary for the 
contaminant to travel a distance L = 50 cm is T m 14.2 hours — probably 
too long to be practical. A suggestion is to reduce the channel length or 
increase the longitudinal flow speed, while increasing the cross-flow velocity 
to be closer to the optimum. 



A Computer algebra handles all the details 



15 




Figure 3: Function R(V) characterising effectiveness of separation of two 
different contaminants for parameters given in Table [l]. 
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A Computer algebra handles all the details 

Just one of the virtues of this centre manifold approach to modelling is that 
it is systematic. This enables relatively straightforward computer programs 
to be written to find the centre manifold and the evolution thereon [12, e.g.]. 

For this problem the iterative algorithm is implemented by a computer 
algebra program written in reduce [] Although there are many details in 
the program, the correctness of the results are only determined by driving to 
zero (line 48) the residual of the governing differential equation, evaluated on 
line 41, to the error specified on line 38 and with boundary and amplitude 
conditions checked on lines 50-52. The other details only affect the rate of 
convergence to the ultimate answer. 



*At the time of writing, information about reduce was available from Anthony 
C. Hearn, RAND, Santa Monica, CA 90407-2138, USA. E-mail: reduce@rand.org 
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1 COMMENT Use iteration to form the centre manifold model of shear 

2 dispersion in a channel with a constant cross-flow of velocity -v. 

3 Flow between y=0 and y=V where V=(b v)/k, rsig=k/nu, A=(v/u*)~2, 

4 u*=-l/2 dp/dx (b k)/(rho mi v) . The centre manifold is parameterized 

5 with c(x,0,t) such that the corrections satisfy c'(x,0,t)=0, 

6 dc' (x,0,t)/dy=0 ; 
7 

8 7» formating for printed output 

9 on div; off allfac; on revpri; factor ev,rsig,df ,a; 

10 7, ev(y) denotes exp(-y) 

11 operator ev; 

12 let { df (ev(y) ,y)=>-ev(y) , ev(0)=>l, ev(V)=>ep>; 

13 ep:=l-l/m; 

14 % operator for solvability, where m=l/(l-exp(-V) ) 

15 operator intv; linear intv; 

16 let { intv(l,y) => V 

17 , intv(y,y) => V~2/2 

18 , intv(y~~q,y) => V" (q+1) / (q+1) 

19 , intv(y~~q*ev(y) ,y)=>-V~q*ep+q*intv(y~ (q-1) *ev(y) ,y) 

20 , intv(y*ev(y) ,y)=> 1-(1+V)*ep 

21 , intv(ev(y) ,y)=> 1/m }; 

22 7, operator to solve d~2h/dy~2+dh/dy = rhs 

23 operator linv; linear linv; 

24 let{linv(l,y) => y-l+ev(y), 

25 linv(y,y) => -y+y~2/2+l-ev(y) , 

26 linv(y-~q,y)=>y-(q+l)/(q+l)-q*linv(y~(q-l) ,y) , 

27 linv(ev(y) ,y)=>-(y+l)*ev(y)+l, 

28 linv(y*ev(y) ,y)=>-(y~2/2+y+l)*ev(y)+l, 

29 linv(y~~q*ev(y) ,y)=>-y~(q+l)/(q+l)*ev(y)+q*linv(y~(q-l)*ev(y) ,y)}; 

30 7o linear solution and velocity profile 

31 depend c,x,t; 

32 let df(c,t)=>g; 

33 7. u:=-2/rsig/V*(y-V*(l-exp(-y*rsig))/(l-exp(-V*rsig))) ; 

34 u : =y* (y-V) * (2*y-6/rsig-V) *rsig/6/V; 

35 h:=c*ev(y) ; 

36 g:=0; 

37 7. iteration, for small d/dx and small reciprocal Prandtl number 

38 let {df (c,x,~q)=>0 when q>2, rsig~2=>0}; 

39 m:=l; '/, optionally neglect e~(-V) terms 

40 repeat begin 

41 eqn:=df (h,t)+u*df (h,x)-df (h,y)-df (h,y,2)-df (h,x,2)*K~2; 

42 7. solvability 

43 gd:=-intv(eqn,y)*m; 

44 g:=g+gd; 

45 7o concentration field 

46 h:=h+linv(eqn+gd*ev(y) ,y) ; 

47 showtime; 

48 end until (eqn=0) ; 

49 7d confirm boundary and amplitude conditions 
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50 eqn:=sub(y=0,h+df (h,y)) ; 

51 eqn:=sub(y=V,h+df (h,y)) ; 

52 eqn:=sub(y=0,h)-c; 

53 % output to file 

54 off nat ; on list; out "sffp.out"; 

55 cmean:=intv(h,y)/V; g:=g; h:=h; 

56 shut "sffp.out"; on nat; 

57 showtime; 

58 end; 



B Weak and moderate cross-flows 

For completeness we record here the model for the case of relatively slow 
cross-flow, or equivalently of relatively high diffusivity. This provides results 
for all parameters V, not just the large values described earlier. 

The non-dimensionalisation used in the main body of this paper is inap- 
propriate in the case of small cross-flow rates. At higher rates the contam- 
inant is restricted to the boundary layer, but in low cross-flow it is spread 
over the channel height. Thus in the case of weak cross- flows we adopt the 
following scalings typical of those used for shear dispersion |8|, [16], e.g.]: 

V — T > * = T? , x = —To, u =-, v = — , p = . (47) 

Quantities are scaled: y with the channel width; t with a cross-channel dif- 
fusion time, r = b 2 /n ~ 1.25 x 10 4 sec; x with the downstream advection 
distance in a cross-channel diffusion time, £ = ur m 1.25 x 10 3 cm; u with 
the mean downstream velocity; and v with a cross-stream diffusion speed, 
n/b ~ 4 x 10~ 4 cm/s. As before V = v^b/n is the main parameter and is 
used to denote a non-dimensional cross-flow velocity, though it may well be 
thought of as an effective channel width, or as the inverse of the molecular 
diffusivity. 

Then after substituting (|4"T|) into the Navier-Stokes and continuity equa- 
tions and dropping stars the equation for the steady horizontal ve- 
locity component u(y) becomes 

V du d 2 u , , , 

— = 12 + — , s.t. u = at y = and y = 1 (48) 

a ay ay 1 

with the nearly parabolic solution 

12a (l- e- Vy/a \ 



u(y) 



V \ 1 - e- v /° y 
y(l-y)U + (l-2y)-+0[^r) I . (19) 
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The advection-diffusion equation for the contaminant becomes 

dc dc ydc d 2 c 1 d 2 c 
dt dx dy dy 2 U 2 dx 2 ' 

where U = ubj k is a downstream Peclet number as before, and with boundary 
conditions 

dc 

Vc + — = at y = and y = 1. (51) 
dy 

In the absence of any x-variations the steady solution is 

Co = Ce~ Vy , (52) 

where as before C = c(x, 0, t) is the concentration of the contaminant at the 
lower plate. The other x-independent solutions, all decaying, are 

c n = C n [V sin (nny) — 2nn cos (niry)} exp (—-^- + A n ^ , (53) 



2 



for n = 1, 2, . . . where the decay rate is 



X n = -^-n 2 n 2 . (54) 

For small V the decay is dominated by the second term above due to cross- 
channel diffusion. The slowest rate of decay to the centre manifold comes 
from the n = 1 mode. Using arguments similar to those given in Section |3] we 
deduce that in the case of small cross-flow rates C\ oc 1 and, consequently, 
the dimensional decay time is expected to be t/|Ai| ~ t/tt 2 w 20min which 
is an order of magnitude larger than that for the strong cross-flows considered 
earlier. This reaffirms the existence of an attractive centre manifold for slowly 
varying solutions, albeit attractive on a larger time scale. 

As before, an iterative procedure was implemented in computer (not 
listed) to solve the contaminant transport equations fl50l)-(|5l]) by finding 
the centre manifold and the evolution thereon (22). The resulting expression 
for the concentration field is 



dx V 2 



( 6 _3 V + 2Vy)(l-£)-^(l-^ 
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The general expressions for the coefficients of the evolution equation (|33| ) are 
quite involved and for brevity here we neglect terms inversely proportional 
to the Prandtl number cr since it is typically small (er ~ 5 x 10 4 for example): 

U = ^(2m-l)-- + 0(a- 1 ) , (56) 

1 e~ v e~ v 4m 2 e~ v — 1 

D = — + 24m— - 144m 2 (2m- 1)— -72 



W V 2 v ' V 3 V 4 

2m- 1 2016 / i \ 

- 720 n^ + - yr + ° ) ' (57) 

e -v 439 , N e -v 

£ = -24m 2 (2m- 1)— + m 2 ( 20m 2 e~ v + 3 ) — 

e" 2V 4 e- 2V 



- 6912m 4 (2m - 1) — -- 24192m 

o«,/o .^SmV^ + S K10 ,48m 2 e- V - 17 

- 864(2m - 1) — 5184 — 

2m - 1 1672704 / ,\ 

- 642816 n^ + + ° (* ) ' (5f 

16m 2 (6mV v + l) %— - i^m 2 (2m - 1) (l5mV v + l)^f 

288 6 ^ 

+ — m 2 (2m - 1) (l68m 2 (100m 2 e~ v + 21)e~ v + 37) — 

„-v 

- 1152m 2 (2m - 1) (24m 2 (l5m 2 e" v - l) e~ v 



V 7 



6912 Q /„„ o / — o _u ~\ _u _„\ a 
m 



.... (60m 2 (25mV v - 2) e~ v - 79 y 
- 6912m 2 (2m - 1) (630mV v + 17 



V 9 

m 2 f 1806m 2 e~ v + 197) e~ v - 33 

- 13824^ — 1 

. 106m 2 e~ v + 29 463m 2 e" v - 237 

- 518400(2m - 1) — - — 829440 — 

V V 

2m - 1 2947995648 , n / A 

- 1208742912^^ + — + O (a' 1 ) . (59) 



These coefficients are plotted in Figure f|. All of them eventually decrease 
in magnitude with increasing cross-flow. The maximum of the effective dif- 
fusion coefficient D is reached at V ~ 3. Thus the "zone broadening" |7| 
associated with the dispersion of the contaminant affects the distribution of 
the contaminant at a greater degree when the cross-flow is relatively weak. 
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Figure 4: The coefficients of the evolution equation fl33|) as functions of the 
cross-channel Peclet number V for the infinite Prandtl number a. 
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The fourth order coefficient F is negative for V < 5 and consequently equa- 
tion (0), in contrast to the case of strong cross- flows reported in Section |], 
predicts stable (decaying) in time evolution of the average concentration of 
the contaminant for all longitudinal wavenumbers. 

The small V expansions of expressions (|55|) -(|59|) are 



U 
D 
E 
F 



C 



1-Vy + -V 2 y 2 
Vy 



6 



dC_ 2 

dx ^ 



+ ^ 12%' - iby + 20 + - [6y' - 15y 



10 



+ 
+ 



120 

V 3 y 



20y A - 36y 3 + 15y 2 - 1 



1 

a 



10y 4 - 24y 3 + 15y 2 + 1 



1260 



60y 4 - 105y 3 + 42y 2 - 7 + 



1 (45^ 



105y 3 + 63y 2 + 7 



+ o 



d 2 C 
dx 2 



^-2 ->)4 



60 2520 



V 6 



100800 



1 1 
W + 210 



1 



60 



V 2 



89 
7920 



V 4 + 



239 
386100 



69300 



1 + l^v 2 - —V 4 + 71629 y 



1365 



2252250 V 
+ 0(V 1 ,V 8 



26879 
85680 



910 
V 2 



2570400 



26311969,^4 1187149277 ^ f 



68372640 



15041980800 



(60) 

(61) 
(62) 
(63) 

(64) 



The expansions for large cross-flow V are 
c = Ce~ Vy 



U 
D 
E 
F 



i <9C 2 _ Vy 

+ IS - ^ e 

ax 



'd 2 C 
ax 2 ' 
2 



2,(7 



2y 



3 6 
V " + V 2 



in 
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a , e 
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4320 

456192 
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1936\ 
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71084 



V 
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The above expressions are equivalent to, but appear a little different from, 



the leading terms in expressions (p0|)- (|32|) , ([34]) and ( p5|) because of the 
different non-dimensionalisation. These large V expressions evidently give 
the behaviour of the coefficients for non-dimensional cross-flow V bigger than 
about 5-10. 
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